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The statistical mechanics of flexible two-dimensional surfaces (membranes) appears in a wide variety of physical 
settings. In this talk we discuss the simplest case of fixed-connectivity surfaces. We first review the current 
theoretical understanding of the remarkable flat phase of such membranes. We then summarize the results of a 
recent large scale Monte Carlo simulation of the simplest conceivable discrete realization of this system jj]]. We 
verify the existence of long-range order, determine the associated critical exponents of the flat phase and compare 
the results to the predictions of various theoretical models. 



1. INTRODUCTION 

Physical membranes, or 2-dimensional surfaces 
embedded in R 3 , are believed to have a high- 
temperature crumpled phase and a low tempera- 
ture flat phase S . The flat phase is characterized 
by long range orientational order in the surface 
normals. Since long-range order is highly unusual 
in 2-dimensional systems, it is worthwhile devel- 
oping a thorough understanding of this phase. 

There are several experimental realizations of 
crystalline surfaces. Inorganic examples are thin 
sheets (< 100 A) of graphite oxide (GO) in an 
aqueous suspension and the rag-like struc- 
tures found in M0S2 

There are also remarkable biological examples 
of crystalline surfaces such as the spectrin skele- 
ton of red blood cell membranes. This is a 
two-dimensional triangulated network of roughly 
70,000 plaquettes. Actin oligomers form nodes 
and spectrin tetramers form links || . Crystalline 
surfaces can also be synthesized in the laboratory 
by polymerising amphiphillic mono- or bi-layers. 
For recent reviews see . 

In this contribution we first review briefly the 
current analytical understanding of rigid mem- 
branes and then summarize the results of a recent 
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large-scale Monte- Carlo simulation 0. 

A Landau-Ginzburg- Wilson effective Hamilto- 
nian for a D-dimensional elastic manifold with 
bending rigidity, embedded in <i-dimensional 



space is 
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where the order field <fi a = d a r, and f(a) is a vec- 
tor in R d . We are clearly dealing with a matrix 
c/) 4 theory. In the crumpled phase (t > 0) the 
position vector r scales with system size L like 
f ~ L v , with v < 1, and hence <p a ~ IS^ 1 . The 
exponent v — 2/du is the size or Flory exponent, 
where is the Hausdorff dimension. An expan- 
sion in <j) a and its gradients is therefore justified. 
In the flat phase (t < 0), the Hamiltonian is sta- 
bilized by the anharmonic terms. In mean field 
theory one finds that the induced metric g a p has 
non-zero expectation value, {d a f ■ dpr) oc b a p- 

The Hamiltonian Eq. (|l|) takes the following 
form in the flat phase 
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where the strain tensor u a p = \{d a f ■ dpr — & a p). 
The upper critical dimension for this model is 4 




r: position 
s: rest position 
u: strain vector 
h: height 




Figure 2. A snapshot of the surface in the flat 
phase for L = 46. 



Figure 1. The Monge gauge. 



and it may be analyzed in a e = 4 — D expan- 
sion. The most pressing physical question is the 
value of the lower critical dimension Di. In par- 
ticular one would like to know if Di < 2, in which 
case physical membranes would indeed have a sta- 
ble flat phase. The Hamiltonian Eq. (||) and the 
strain tensor u a f3 can be re-written in the Monge 
gauge (see Figure |l|) as 



H *ff = \J d D a[ K (d 2 hf+2^ 
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and 



u aP = -^(daUp + dpu a + d a hdph), 



(3) 
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where /i and A are the Lame coefficients.^ 

By rescaling a — > sa, one sees that Eq. (||) in 4 
dimensions is a function only of the dimensionless 
parameters 

U = -^r and X = —. (5) 

K K 

Aronovitz and Lubensky (AL) determined the 
RG flow of ft and A within the e-expansion, at 
fixed co-dimension, and found a globally attrac- 
tive IR-stable fixed point at infinite bending rigid- 
ity [[l0| . This fixed point should control the prop- 
erties of the whole flat phase pl|-p^[. In par- 
ticular, we can introduce the anomalous scaling 

2 Note that it is impossible for the surface to trade stretch- 
ing energy for bending energy so as to make u a p every- 
where zero, since the two phonon degrees of freedom are 
insufficient to cancel the three independent components of 
the strain tensor. 



dimensions for the running coupling constants, 
Kfl(g) ~ q~ v and ^r{o) ~ A fl (g) ~ q Vu - The 
generalization of Eq. (|^) to D dimensions relates 
the exponents i] and rj u with the scaling identity 



i lu = 4 - D - 2rj. 



(6) 



The roughness exponent £, which governs the 
scaling of the height-height correlation function, 
is related to the exponent r\ by the scaling rela- 
tion £ — (4 — D — r))l2. One may also study the 
Hamiltonian of Eq. (|3|) in a large-d expansion in 
the non-linear sigma model limit (infinite Lame 
coefficients) [Q. Finally one can solve a set of 
self-consistent equations (SCSA) for the scaling 
exponents of the renormalized coupling constants 
k, /i and A The predictions for (, ry, r\ u and 
v from the different approximations, and the re- 
sults of our Monte Carlo (MC) simulations are 
summarized in Table [l| 

2. THE MODEL 

In the simplest discretized version the crys- 
talline surface is modeled by a regular triangu- 
lar lattice with fixed connectivity, embedded in d- 
dimensional space. Typically the link-lengths of 
the lattice are allowed some limited fluctuations. 
This is usually modeled by tethers between hard 
spheres or by introducing some confining pair po- 
tential with short-range repulsion between nodes 
(monomers), such as Lennard- Jones. In some 
cases, the bending energy is explicitly introduced, 
and is represented by a ferromagnetic-like inter- 
action between the normals to nearest-neighbor 



3 



Table 1 

Theoretical predictions and numerical results. 
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Vu 


v = 2/d„ 
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AL 


13/25 or 0.52 


2/25 or 0.08 


1 


24/25 or 0.96 


Large-d 


2/3 


2/3 


1 


2/3 


SCSA 


0.590 


0.358 


1 


0.821 


MC 


0.64(2) 


0.50(1) 


0.95(5) 


0.62 



"plaquettes" g. 

We are interested in a much simpler model of a 
crystalline surface, inspired by the Polyakov ac- 
tion for strings |16|. The tethering potential be- 
tween the particles is a simple Gaussian potential, 
with vanishing equilibrium length. Since the equi- 
librium length defines the bare elastic constants 
of the surface one may wonder whether such a 
model indeed exhibits the same phase diagram, 
and, in particular, if it has a stable flat phase. 

N particles are arranged in a regular triangular 
mesh. Each node in the interior has 6 neighbors. 
The action is composed of a tethering potential 
and a bending energy term: 



rib) 



(7) 



(ab) 



where r is the position of node i, and n a is the 
unit normal to face a. The sums extend to near- 
est neighbors. We point out that the normal- 
normal interaction translates to a next-to-nearest 
neighbor interaction (like a V 2 term.) We do 
not include any minimum distance between the 
nodes. The model describes a phantom surface, 
since there is no self-avoidance term. While self- 
avoidance changes the nature of the crumpling 
transition, in the flat phase it is irrelevant (see Y. 
Kantor in ||.) 

We choose to simulate a surface with free 
boundaries, since it simplifies the analysis of the 
correlation functions. The trade off is that we 
have to take careful account of edge fluctuations. 
This is described and illustrated in detail in 

3. OBSERVABLES 

In order to estimate the exponents shown in Ta- 
ble we measured several observables. For finite 
sizes, the specific heat Cy peaks in the vicinity 




<1 + K) 



Figure 3. The specific heat Cy versus k/(1 + k). 



of a second order phase transition, or at coupling 
k — k c . The peak diverges with the exponent 
a ~ 0.4. We are currently increasing the statis- 
tics around the phase transition in order to de- 
termine a with greater accuracy. The main un- 
certainty affecting the estimate for a comes from 
the estimate of k c . We also measure Cy in or- 
der to locate a region of the flat phase suitable 
for studying its scaling behavior. Note that, for a 
surface of finite size, the bending rigidity controls 
the importance of finite size effects. In order to 
obtain reliable finite size scaling one has to tune 
the correlation length £ ~ L. 

While Cy indicates the location of the transi- 
tion, it tells little about the nature of the phases 
on each side. Thus we measured the shape tensor 
and computed its eigenvalues. The shape tensor 
is the off-diagonal part of the inertia tensor I a @ , 
being in some sense orthogonal to it. It is defined 
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Figure 4. The distribution of eigenvalues of the 
shape tensor. 



as 

S af 3 = (^2r a (a)r (ayj (8) 

where a, refer to the components of r, and the 
subscript c indicates connected expectation val- 
ues. The radius of gyration R 2 g = tiS changes 
drastically across the transition. While for k < k c 
(hot phase) R g has small values compared with 
the linear size L of the surface, for k > n c R g 
is large. More importantly the finite size scaling 
of R g ~ L" defines the size exponent v — 2/oIh- 
In the crumpled phase R g ~ log(L), and v = 0. 
At the critical point ^ has a non-trivial value. 
Above the critical point, in the flat phase, v — 1 
or <1h = 2. Our measurement of v = 0.95(5) at 

K = 1.1. 

The eigenvectors of S define a body-fixed frame 
on the surface; the eigenvalues of S are the aver- 
age square dispersions in the direction of the as- 
sociated eigenvalue and tell more about the shape 
of the surface. Figure || shows the distribution of 
eigenvalues p(X) in the three regions of the phase 
diagram. Box a) shows p in the crumpled phase: 
all three eigenvalues have identical distribution 
and the system is isotropic. Box b) shows p in 
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Figure 5. The roughness exponent C versus the 
diameter of the hexagonal sub-set of the surface. 



the vicinity of the transition: the system is still 
isotropic but the eigenvalues have fluctuations on 
large scales. Box c) shows p in the flat phase: the 
surface is no longer isotropic — there is a well de- 
fined thickness and lateral extension, correspond- 
ing, respectively, to the left and right peak. 

The roughness exponent £ is measured from the 
finite size scaling of the average thickness of the 
surface. The minimum eigenvalue of S provides 
just this information. This particular observable 
is very sensitive to boundary effects. In fact if 
an edge of the surface is "curled" , the eigenvalue 
will be considerably larger even if the surface is 
locally quite smooth. 

In order to determine the influence of the 
boundary, we measured the average square thick- 
ness of several concentric hexagonal sub-sets of 
the surface of varying diameter D. We found that 
C plateaus in the range L/4 < D < 3L/4, shows 
boundary effects for D > 3L/4 and discretization 
effects for D < L/4. We quote the value extracted 
from the plateau region (see Figure ||) . 

The exponent rj u can be extracted from the 
phonon fluctuations. It is in fact straight-forward 
to relate this observable and the exponent using 
the finite size scaling relation 

<|u| 2 )~Z7K (9) 
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Table 2 

The number of thermalized sweeps collected per data point in the flat phase. The last column indicates 
the autocorrelation time for the slowest mode in the system, the radius of gyration. 
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K = 1.1 


n = 2.0 




32 


31 x 10 e 


26 x 10 B 


3 x 10 4 


4G 


51 x 10 6 


42 x 10 6 


7 x 10 4 


64 


47 x 10 6 


44 x 10 6 


1.2 x 10 5 


128 


74 x 10 6 




5 x 10 5 



We have performed precise measurements of the 
phonon fluctuations, described in detail in ref. 
jlj. The determination of this exponent is im- 
portant since it provides and independent consis- 
tency check of the scaling relation (||). Our mea- 
surement, shown in Table El is in good agreement 
with the theoretical predictions and the scaling 
relation. 

Our final measurement of the properties of the 
flat phase of this model is the normal-normal cor- 
relation function. We expect the correlation func- 
tion to fall off to a non-zero asymptote like 

(n a -n ) ~ C+-1 (10) 

where r is the geodesic distance between the cen- 
ter o and the point a. Since the boundaries are 
free, the correlation function is not translation- 
ally invariant: we therefore fix the origin at the 
center of the surface and discard data too close 
to the boundary. 

In Figure § we show the behavior of the cor- 
relation function for different values of the lat- 
tice size. From the data in the figure it is clear 
that the normal-normal correlation function has a 
non-zero asymptote. The value of the asymptote 
grows with system size, and we believe it has a 
non zero infinite volume limit. The fit of the data 
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to Eq. 10 gives yet another independent estimate 
of the exponent r\. Our best fit to the L = 128 
data gives rj ~ 0.62, in good agreement with our 
previous estimates. 

In conclusion, we have shown that the simple 
gaussian model of Eq. [j] faithfully reproduces the 
expected critical behavior of the AL fixed point. 
The relative simplicity of this model enables us 
to simulate surfaces of realistic size. The spectrin 
network of a blood cell has about 70,000 plaque- 
ttes; our largest surface has 32,768 plaquettes. 
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Figure 6. The normal-normal correlation function 
for different lattice sizes in the flat phase. The 
dashed lines are our best fit to Eq. nfl. 



4. NUMERICAL METHODS 

We will now turn our attention to the technical 
aspects of the simulation methods and the com- 
puters used. We use a Monte Carlo algorithm 
with a local Metropolis update. We choose a new 
position for a given node in a box of size e 3 cen- 
tered on the old position, and we accept it ac- 
cording to the standard Metropolis test. A single 
Monte Carlo sweep consists of an update of all 
the nodes of the surface. The box size e is ad- 
justed so as to keep the acceptance ratio around 
50%. 

Crystalline surfaces are characterized by ex- 
tremely long autocorrelation times. The modes 
which suffer most are the ones related to the 
global shape of the surface in the embedding 
space. In table |^ we show the amount of data 
collected at various sizes/bending rigidities, and 
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Table 3 



Wall-clock timing results of the benchmark runs for different serial and parallel codes. The number are 
in seconds. 





Serial code 




MPI code 


HPF code 
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B 


A 
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Serial 


334 






4800 




2 




705 


1410 


2189 


4378 


4 




400 


1600 


2375 


9500 


8 




212 


1696 







the corresponding autocorrelation time for the ra- 
dius of gyration. We estimate the critical slowing 
down exponent z to be ~ 2, as expected for a 
local algorithm. 

We are currently investigating several ways of 
improving the Monte Carlo simulations meth- 
ods. In particular, we are studying different al- 
gorithms, like over-relaxation and uni/multi-grid 
Monte Carlo. Previous numerical studies have 
successfully taken the advantage of Fourier accel- 
eration [p^|-^9|, and it is argued that multi-grid 
methods should perform similarly well |^0|, p. 41]. 
Preliminary tests of these algorithms give very en- 
couraging results: over-relaxation reduces signifi- 
cantly the autocorrelation time, while we expect 
the multi-grid to improve the dynamical exponent 
z. 

We used several workstations for most of the 
runs, but we simulated the largest lattice size 
(L = 128) using a MASPAR MPI, a massively 
parallel processor, with exactly 16,384 nodes. 
Due to its hardware, the MPI is best suited 
to run a Monte Carlo simulation with local up- 
date. Massively parallel processors are no longer 
common: workstation clusters with high speed 
switches are becoming more and more popu- 
lar. This trend is reflected in the parallclization 
strategies that we are currently investigating. 

While High Performance Fortran (HPF) is 
clearly the better choice for massively parallel 
computers, such as the MPI, it is still unclear 
whether its efficiency (or lack thereof) justifies its 
use in clustered environments. The alternative is 
to use the Message Passing Interface (MPI) . This 
is a library of functions and a set of programming 
models which allows one to distribute a simula- 
tion over a cluster of workstations. 



It is much simpler to write programs in HPF 
than MPI, since the HPF compiler automatically 
distributes the data on the various processors 
and it introduces the appropriate parallel instruc- 
tions. The trade-off is that the programmer has 
little control over how the parallelization is actu- 
ally done, and often the compiler does not exploit 
all of the potential parallelism of the algorithm. 
MPI allows much more flexibility and tailoring, 
but it is much harder to use since the paralleliza- 
tion must be coded by hand. 

We have performed a series of benchmarks, 
using a simplified version of the code, in order 
to establish the efficiency of the different par- 
allclization techniques. In Table || we compare 
the results of the benchmark runs. The first col- 
umn shows the number of processors used in the 
run. This is an important factor: more proces- 
sors mean more computing power, but also more 
message passing between the processors. While 
massively parallel computers had relatively low 
latency times for communication, clustered work- 
stations often rely on a common switchboard with 
high latency times. The first column lists the 
elapsed time for a single processor run of the 
scalar code — this is the reference. The third 
and fourth column have the timings for the MPI 
code, written in c. Column A shows the elapsed 
time, while column B shows the integrated time 
(wall-clock time x number of processors). The 
fifth and sixth column show the timing of the 
HPF code. We note that, while for the MPI code 
the integrated time (column 3) grows slowly with 
the number of processors, for the HPF code (col- 
umn 6) the time doubles, indicating very ineffi- 
cient message passing. Based on our limited ex- 
perience with our particular code, we believe that 
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HPF is still far from being a viable choice for par- 
allel programming on clustered environments. 

Finally, we believe that farming (or running in- 
dependent serial simulations on many processors) 
is still the most efficient solution for simulations 
small enough to fit on a single processor, since the 
integrated time for the MPI code is still higher 
than the single processor time. 
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